library(data.table)
library(magrittr)
library(parallel)
library(dplyr)
library(boot)
library(INLA)
library(progress)
sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49")
sel1=paste(sel1,".snp",sep="")
snpdiff=c()
pb <- progress_bar$new(total = 100)
for (i in 1:length(sel1)){
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="T",]
result=matrix(,nrow(file),ncol=2)
if(dim(file)[1]!=0){
fn="normal_tumor"
for(k in 1:2){
files=file[,c(paste(unlist(strsplit(fn,"_"))[k],"_reads1",sep=""),paste(unlist(strsplit(fn,"_"))[k],"_reads2",sep=""))]
for(j in 1:nrow(files)){
reads=as.numeric(as.matrix(files)[j,])
ref=reads[1:(ncol(files)/2)]
var=reads[(ncol(files)/2+1):ncol(files)]
df=data.frame(y=var,Ntrials=ref+var)
null <- logit(median(df[,1]/df[,2]))
formula = y ~ 1
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
m <- m1$marginals.fixed[[1]]
	lower_p <- inla.pmarginal(null, m)
	upper_p <- 1 - inla.pmarginal(null, m)
	post_pred_p <- 2 * (min(lower_p, upper_p))
	coef <- m1$summary.fixed
	result[j,k]=post_pred_p
	}
}
print(c(paste(i,"is done")),sep="")
file$normal_pvalue=result[,1]
file$tumor_pvalue=result[,2]
sels=which(file$normal_pvalue<=0.05 |file$tumor_pvalue<=0.05)
fila=file[sels,]
countdiff=paste(fila$chrom,fila$position,fila$ref,fila$var,sep="_")
snpdiff=unique(c(snpdiff,countdiff))
}
  pb$tick()
  Sys.sleep(1 / 100)
}


###以下为合并文件
rt=data.frame(as.matrix(snpdiff,1,length(snpdiff)))
names(rt)=c("unitID")
rt=tidyr::separate(rt,unitID,into=c("chr","pos","ref","var"),sep="_")
rt$unitID=paste(rt$chr,rt$pos,sep=":")
positions=as.character(unlist(rt$unitID))

i=1
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="T",]
if(dim(file)[1]>0){
file$unitID=paste(file[,1],":",file[,2],sep="")
rownames(file)=file$unitID
final=file[positions,c("normal_reads1","normal_reads2","normal_var_freq","tumor_reads1","tumor_reads2","tumor_var_freq","unitID")]
final=final[complete.cases(final$unitID),]
final$tumor_var_freq=as.numeric(gsub("%","",final$tumor_var_freq))/100
final$normal_var_freq=as.numeric(gsub("%","",final$normal_var_freq))/100
final$total_reads=as.numeric(final$normal_reads1)+as.numeric(final$normal_reads2)+as.numeric(final$tumor_reads1)+as.numeric(final$tumor_reads2)
sels=which(final$total_reads>=10)
final=final[sels,]
final=final[,-8]

colnames(final)=c(paste(strsplit(gsub(".snp","",sel1[i]),"_")[[1]][1],c("_reads1","_reads2","_var_freq"),sep=""),paste(strsplit(gsub(".snp","",sel1[i]),"_")[[1]][2],c("_reads1","_reads2","_var_freq"),sep=""),"unitID")
}

for (i in 2:length(sel1)){
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="T",]
if(dim(file)[1]>0){
file$unitID=paste(file[,1],":",file[,2],sep="")
rownames(file)=file$unitID
final1=file[positions,c("normal_reads1","normal_reads2","normal_var_freq","tumor_reads1","tumor_reads2","tumor_var_freq","unitID")]
final1=final1[complete.cases(final1$unitID),]
final1$tumor_var_freq=as.numeric(gsub("%","",final1$tumor_var_freq))/100
final1$normal_var_freq=as.numeric(gsub("%","",final1$normal_var_freq))/100
final1$total_reads=as.numeric(final1$normal_reads1)+as.numeric(final1$normal_reads2)+as.numeric(final1$tumor_reads1)+as.numeric(final1$tumor_reads2)
sels=which(final1$total_reads>=10)
final1=final1[sels,]
final1=final1[,-8]
colnames(final1)=c(paste(strsplit(gsub(".snp","",sel1[i]),"_")[[1]][1],c("_reads1","_reads2","_var_freq"),sep=""),paste(strsplit(gsub(".snp","",sel1[i]),"_")[[1]][2],c("_reads1","_reads2","_var_freq"),sep=""),"unitID")

final=merge(final,final1,by="unitID",all=T)
}
}

###########读取文件并进行INLA分析

######ASE function要求结果file文件为count数，且顺序为
######twin1_control_ref,twin1_case_ref,twin2_control_ref,twin2_case_ref,
######twin1_control_var,twin1_case_var,twin2_control_var,twin2_case_var.

sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49")
#sel1=c("X2B_X1T","M8_M7","M50_M49","M52_M51");sel2=c("M28_M27")
#sel2=c("M30_M29","M26_M25","M35_M36")
#sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M42_M41","M44_M43","M48_M47","M50_M49","M52_M51","M12_M11");sel2=c("M28_M27","M30_M29","M26_M25","M36_M35")  ####包括抑郁症
sel3=c("M18_M17","M20_M19","M22_M21","M40_M39")

filekp=read.table("../bayes/AShM.txt",head=T,sep = "\t")
ref=paste(unlist(strsplit(sel1,"_")),"_reads1",sep="")
var=paste(unlist(strsplit(sel1,"_")),"_reads2",sep="")
file=filekp[,c(ref,var)]

ASE=function(x)
{
if (length(which(!is.na(x)))/4>=1){
sel=which(!is.na(x))
rct=x[sel]
ref=rct[1:(length(rct)/2)]
var=rct[(length(rct)/2+1):(length(rct))]
if (length(sel)/4>1){
df=data.frame(y=var,Ntrials=ref+var,x1=rep(c(0,1),length(sel)/4),x2=rep(c(0:(length(sel)/4-1)),each=2))
#header=c(rep("y",length(sel)/2),rep("Ntrials",length(sel)/2),"x1","x2")
#colnames(df)=header
formula = y ~ 1 + f(x2, model = "iid") + x1
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
formula = y ~ 1 + f(x2, model = "iid")
m0 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
b10 <- exp(m1$mlik[2] - m0$mlik[2])
}else{
df=data.frame(y=var,Ntrials=ref+var,x1=rep(c(0,1),length(sel)/4))
#header=c(rep("y",length(sel)/2),rep("Ntrials",length(sel)/2),"x1")
#colnames(df)=header
formula = y ~ 1 + x1
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
formula = y ~ 1
m0 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
b10 <- exp(m1$mlik[2] - m0$mlik[2])
}
return(b10)
}
}

Sys.time()
rt=apply(file,1,ASE)
Sys.time()
result=as.numeric(as.character(rt))

filekp$Bayes_in_DC=result

